Многомерные модели авторегрессии (VAR)

На практике нас чаще всего интересуют процессы многомерные, поведение которых в некотором смысле связано между собой. Итак пусть теперь \(x_{t}\)-векторный случайный процесс размерности \((n×1)\). Векторным процессом авторегрессии порядка \(p\) (сокращенно VAR(p) ) будем называть процесс вида \[x_{t}=с+\Phi_{1}x_{t-1}+...+\Phi_{p}x_{t-p}+\epsilon_{t}\]

где c- вектор (n×1) констант, \(\Phi_{i}; i=1,...p\)- матрицы размера (n×n) и \(\epsilon_{t}\)-многомерный (n×1)процесс белого шума

\[E\epsilon_{t}=0\] \[E\epsilon_{t}\epsilon_{s}=\Omega\] при \(t=s\); 0 иначе, где \(\Omega\)- симметричная положительно определенная (n×n) матрица ковариаций.

Выпишем соотношение, например, для первого из компонент вектора \(x_{t}\) \[x_{1,t} = c_1+\Phi_{1,1}^{(1)}x_{1,t-1}+\Phi_{1,2}^{(1)}x_{2,t-1}...+\Phi_{1,n}^{(1)}x_{n,t-1}+ \Phi_{1,1}^{(2)} x_{1,t-2}+\Phi_{1,2}^{(2)} x_{2,t-2}+...+\Phi_{1,n}^{(2)} x_{n,t-2} +...+\] \[\Phi_{1,1}^{(p)} x_{1,t-p}+\Phi_{1,2}^{(p)}x_{2,t-p}...+\Phi_{1,n}^{(p)}x_{n,t-p}+\epsilon_{t,1}\]

Таким образом каждая из компонент процесса представляет собой одновременно как линейную регрессионную модель от других компонент процесса, но взятых с задержкой по времени и авторегрессионную модель от самой себя.

Используя векторный оператор задержки \(Bx_{t}=x_{t-1}\) модель может быть записана в форме

\[[I_{n}-\Phi_1B-...-\Phi_{p}B^{p}]x_{t}=c+\epsilon_{t}\]

или

\[\Phi(B)x_{t}=c+\epsilon_{t}\]

где \(\Phi(B)\)- матрица размера (n×n) полиномов от оператора задержки, каждый элемент которой представляет собой скалярный полином от задержки

\[\Phi(B)=[\delta_{ij}-\Phi_{ij}B-...-\Phi_{ij}^{p}B^{p}]; i,j=1,...,n\]

где \(\delta_{ij}=1\) при \(i=j\), иначе 0

Процесс VAR(p) называют стационарным в широком смысле (или просто стационарным), если его первый и второй момент \(Ex_{t}\) и \(Ex_{t}x_{t-j}′\) соответсвенно не зависит от времени \(t\).

Если процесс стационарен, то мы можем взять математическое ожидание от обеих частей и вычислить математическое ожидание процесса. Обозначим вектор математических ожиданий через \(\mu=Ex_{t}\), тогда \[\mu=с+\Phi₁\mu+...+\Phi_{p}\mu\] Откуда \[\mu=(I_{n}-\Phi_1-...-\Phi_{p})^{-1}c\] Используя вектор математических ожиданий \(\mu\), модель VAR(p) может быть эквивалентно переписана в терминах отклонений от среднего \(\mu\). \[x_{t}-\mu=\Phi_1(x_{t-1}-\mu)+...+\Phi_{p}(x_{t-p}-\mu)+\epsilon_{t}\]

Предложение. Многомерная модель VAR(p), так и одномерная модель AR(p) может быть переписана как многомерная модель VAR(1) Доказательство. Действительно, обозначим через \(\xi_{t}\)- многомерный размера (np×1)процесс вида \[\begin{equation*}\xi_{t}=\left[ \begin{matrix}x_t-\mu\\ x_{t-1}-\mu\\ ....\\x_{t-p}-\mu\end{matrix}\right] \end{equation*}\] Пусть \(F\) матрица размера (np×np) \[\begin{equation*}F=\left[ \begin{matrix}\ \Phi_1&\Phi_2&\Phi_3&...&\Phi_p\\ I_n&0&0&...&0\\ 0&I_n&0&...&0\\ ...\\ 0&0&...&I_n&0&\\ \end{matrix}\right] \end{equation*}\] и через \(ν_{t}\)- (np×1)-й процесс \[\begin{equation*}v_{t}=\left[ \begin{matrix}\epsilon_t\\0\\...\\0\end{matrix}\right] \end{equation*}\] Тогда многомерный VAR(p) процесс может быть эквивалентно переписан, как VAR(1) процесс вида \[\xi_{t}=F\xi_{t-1}+ν_{t}\] Доказанный факт чаще всего используется при доказательстве различных утверждений относительно процесса VAR(p).

Условия стационарности процесса VAR(p)

Для процесса VAR(1) \(\xi_{t}\) , введенного выше справедливо следующее выражение\[\xi_{t+s}=ν_{t+s}+Fν_{t+s-1}+F^2ν_{t+s-2}+...+F^{s-1}ν_{t+1}+F^{s}\xi_{t}\]

Для того, чтобы процесс был стационарным влияние шума \(\epsilon_{t}\) в конце концов должно исчезнуть. Оказывается, что если все собственные значения матрицы \(F\) лежат внутри единичного круга, то процесс \(\xi_{t}\) будет стационарным в широком смысле. Без доказательства примем на веру следующее утверждение Теорема. Все \(\lambda\)- собственные значения матрицы \(F\), удовлетворяют уравнению \[|I_{n}\lambda^{p}-\Phi_1\lambda^{p-1}-\Phi_2\lambda^{p-2}...-\Phi_{p}|=0\] Следовательно, процесс VAR(p) является стационарным тогда и только тогда, когда все собственные значения матрицы F по модулю меньше единицы \(|\lambda|<1\), что в свою очередь эквивалентно тому, все корни \(z\) уравнения \[|I_{n}-\Phi_1z-\Phi_2z²...-\Phi_{p}z^{p}|=0\] лежат вне единичного круга. Доказательство этого утверждения можно найти в J.Hamilton. Time Series Analysis p.285.

Автоковариационная функция векторного процесса

Для стационарного векторного процесса \(x_{t}\) размера (n×1) \(k\)-ой автоковариацией назовем матрицу (n×n)\[\Gamma_{k}=E[(x_{t}-\mu)(x_{t-k}-\mu)′]\] Вспомним, что для одномерных процессов автоковариационная функция является четной функцией \(c_{k}=c_{-k}\) , для векторного процесса это не верно. \[\Gamma_{k}\ne\Gamma_{-k}\] Однако, нетрудно убедиться, что \[\Gamma_{k}′=\Gamma_{-k}\]

Действительно. Для стационарного процесса будет справедливо \[\Gamma_{k}=E[(x_{t+k}-\mu)(x_{(t+k)-k}-\mu)′]=E[(x_{t+k}-\mu)(x_{t}-\mu)′]\] Транспонируем \[\Gamma_{k}′=E[(x_{t}-\mu)'(x_{t+k}-\mu)]=\Gamma_{-k}\] В языке R соответствующая библиотека для работы с многомерными авторегрессионными моделями называется mAr Cмоделируем VAR(2) процесс размерности 3

Сначала зададим модель вектор констант \(c\) и матрицы \(\Phi_1,\Phi_2\) и ковариационную матрицу \(\Omega\)

library(mAr)
## Loading required package: MASS
w <- c(1,4,8)
Phi1 <- c(0.1,-0.2,-0.3)
Phi1 <- rbind(Phi1,c(0.1,0.-0.2,0.3))
Phi1 <- rbind(Phi1,c(0.2,0.2,0.2))
Phi1
##      [,1] [,2] [,3]
## Phi1  0.1 -0.2 -0.3
##       0.1 -0.2  0.3
##       0.2  0.2  0.2
Phi2 <- c(0.1,0.1,0.1)
Phi2 <- rbind(Phi2,c(0.15,0.1,0.15))
Phi2 <- rbind(Phi2,c(0.1,0.1,0.1))
Phi2
##      [,1] [,2] [,3]
## Phi2 0.10  0.1 0.10
##      0.15  0.1 0.15
##      0.10  0.1 0.10
Phi <- cbind(Phi1,Phi2)
Phi
##      [,1] [,2] [,3] [,4] [,5] [,6]
## Phi1  0.1 -0.2 -0.3 0.10  0.1 0.10
##       0.1 -0.2  0.3 0.15  0.1 0.15
##       0.2  0.2  0.2 0.10  0.1 0.10
Omega <- diag(4,3,3)
Omega  
##      [,1] [,2] [,3]
## [1,]    4    0    0
## [2,]    0    4    0
## [3,]    0    0    4

Cформируем из матриц \(\Phi_1,\Phi_2\) матрицу \(F\) для того, чтобы проверить стационарность модели

I <- diag(1,3,3)
Zero <- diag(0,3,3)
IZ <- cbind(I,Zero)
F <- rbind(Phi,IZ)
F
##      [,1] [,2] [,3] [,4] [,5] [,6]
## Phi1  0.1 -0.2 -0.3 0.10  0.1 0.10
##       0.1 -0.2  0.3 0.15  0.1 0.15
##       0.2  0.2  0.2 0.10  0.1 0.10
##       1.0  0.0  0.0 0.00  0.0 0.00
##       0.0  1.0  0.0 0.00  0.0 0.00
##       0.0  0.0  1.0 0.00  0.0 0.00
eval <- eigen(F,only.values= TRUE)
modev <- Mod(eval$values)
barplot(modev,col = "blue",main = " Abs value of eigen values")
abline(h=1,col = "black")

Все собственные значения по модулю <1, следовательно ряд стационарен Моделируем и смотрим результат

VAR2 <- mAr.sim(w, Phi, Omega, 100)
matplot(VAR2,type="l",lty = 1, col=c("magenta","blue","green"),lwd=2,main = "Simulated Stationary VAR(2) Model")

Оценим по рядам параметры AR(2) и сравним с теми, что были заданы при моделировании

model <- mAr.est(VAR2,2)
model$wHat
## [1] -0.6347867  3.5109273  6.6697915
model$AHat
##            [,1]       [,2]       [,3]        [,4]        [,5]       [,6]
## [1,] 0.04372861 -0.2509743 -0.1325047 -0.03993169 -0.03283639 0.07787981
## [2,] 0.08610531 -0.3168418  0.3799488  0.07151865  0.05832944 0.17802060
## [3,] 0.26323894  0.2016552  0.2851038  0.18897199  0.12260496 0.13475527
model$CHat
##              [,1]      [,2]         [,3]
## [1,]  6.352889561 0.5562586 -0.007961853
## [2,]  0.556258618 3.8788206  0.661925614
## [3,] -0.007961853 0.6619256  3.278477323

К сожалению метод не дает стандартных ошибок оценок и р-значений проверок гипотез о равенстве нулю параметров модели.

Другой метод оценки находится в библитетеке vars. Здесь уже есть выдача с проверкой гипотез о значимости оценок

В этой библиотеке

library(vars)
## Warning: package 'vars' was built under R version 3.2.5
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 3.2.5
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 3.2.5
## Loading required package: urca
## Warning: package 'urca' was built under R version 3.2.4
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 3.2.4
varsimest <− VAR( VAR2 , p = 2 ,season = NULL, exogen = NULL)
summary(varsimest)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: X1, X2, X3 
## Deterministic variables: const 
## Sample size: 98 
## Log Likelihood: -619.11 
## Roots of the characteristic polynomial:
## 0.5978 0.4694 0.2326 0.2326 0.2032 0.1316
## Call:
## VAR(y = VAR2, p = 2, exogen = NULL)
## 
## 
## Estimation results for equation X1: 
## =================================== 
## X1 = X1.l1 + X2.l1 + X3.l1 + X1.l2 + X2.l2 + X3.l2 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)  
## X1.l1  0.04373    0.10522   0.416   0.6787  
## X2.l1 -0.25097    0.13729  -1.828   0.0708 .
## X3.l1 -0.13250    0.14170  -0.935   0.3522  
## X1.l2 -0.03993    0.10929  -0.365   0.7157  
## X2.l2 -0.03284    0.13535  -0.243   0.8089  
## X3.l2  0.07788    0.14596   0.534   0.5949  
## const -0.63479    1.91503  -0.331   0.7410  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 2.52 on 91 degrees of freedom
## Multiple R-Squared: 0.06887, Adjusted R-squared: 0.007481 
## F-statistic: 1.122 on 6 and 91 DF,  p-value: 0.356 
## 
## 
## Estimation results for equation X2: 
## =================================== 
## X2 = X1.l1 + X2.l1 + X3.l1 + X1.l2 + X2.l2 + X3.l2 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## X1.l1  0.08611    0.08222   1.047 0.297731    
## X2.l1 -0.31684    0.10728  -2.953 0.003998 ** 
## X3.l1  0.37995    0.11073   3.431 0.000905 ***
## X1.l2  0.07152    0.08540   0.837 0.404508    
## X2.l2  0.05833    0.10576   0.552 0.582633    
## X3.l2  0.17802    0.11405   1.561 0.122015    
## const  3.51093    1.49637   2.346 0.021132 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 1.969 on 91 degrees of freedom
## Multiple R-Squared: 0.2573,  Adjusted R-squared: 0.2083 
## F-statistic: 5.255 on 6 and 91 DF,  p-value: 0.0001095 
## 
## 
## Estimation results for equation X3: 
## =================================== 
## X3 = X1.l1 + X2.l1 + X3.l1 + X1.l2 + X2.l2 + X3.l2 + const 
## 
##       Estimate Std. Error t value Pr(>|t|)    
## X1.l1  0.26324    0.07559   3.483 0.000765 ***
## X2.l1  0.20166    0.09863   2.045 0.043780 *  
## X3.l1  0.28510    0.10180   2.801 0.006228 ** 
## X1.l2  0.18897    0.07851   2.407 0.018104 *  
## X2.l2  0.12260    0.09723   1.261 0.210558    
## X3.l2  0.13476    0.10485   1.285 0.201986    
## const  6.66979    1.37571   4.848 5.11e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 1.811 on 91 degrees of freedom
## Multiple R-Squared: 0.3908,  Adjusted R-squared: 0.3506 
## F-statistic: 9.728 on 6 and 91 DF,  p-value: 2.908e-08 
## 
## 
## 
## Covariance matrix of residuals:
##           X1     X2        X3
## X1  6.352890 0.5563 -0.007962
## X2  0.556259 3.8788  0.661926
## X3 -0.007962 0.6619  3.278477
## 
## Correlation matrix of residuals:
##           X1     X2        X3
## X1  1.000000 0.1121 -0.001745
## X2  0.112058 1.0000  0.185619
## X3 -0.001745 0.1856  1.000000

В этой библиотеке также присутствут очень удобный метод VARselect, который позволяет оценить порядок модели \(p\), когда он неизвестен Это делается с помощью 4 информационных критериев, среди которых и известный нам уже критерий Акаике, который в многомерном случае выглядит следубщим образом \[AIC(n) = \ln \det(\tilde{Σ}_u(n)) + \frac{2}{T}n K^2 \quad\] где \[\tilde{Σ}_u (n) = T^{-1} ∑_{t=1}^T \hat{u}_t \hat{u}_t'\] где \(u_t\) вектора остатков.

nfocrit <− VARselect(VAR2, lag.max = 5, type="const")
nfocrit
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      1      1      1      1 
## 
## $criteria
##                1         2          3          4          5
## AIC(n)  4.525521  4.593109   4.716654   4.721484   4.739506
## HQ(n)   4.655874  4.821226   5.042535   5.145130   5.260916
## SC(n)   4.848116  5.157651   5.523141   5.769918   6.029886
## FPE(n) 92.357841 98.880511 112.056382 112.918666 115.494923

Про другие критерии можно прочитать в help по этому методу. Прогнозирование также, как и в ARMA модели это условное математическое ожидание, которое в данном случае при прогнозе на один шаг времени вперед будет

\[E[x(t+1)|x_1,...x_t] = с+\Phi_{1}x_{t}+...+\Phi_{p}x_{t-p+1}\] прогноз на \(l\) шагов вперед вычисляется последовательно через прогнозы на \(l-1,...1\) шагов. Соответствующий метод показан ниже

predictions <− predict ( varsimest , n.ahead = 25 ,ci = 0.95)
matplot(predictions$fcst$X1,type = "b",pch = 21,lty = 1,lwd = 3,main = "Forecast of x1")

matplot(predictions$fcst$X2,type = "b",pch = 21,lty = 1,lwd = 3,main = "Forecast of x2")

matplot(predictions$fcst$X3,type = "b",pch = 21,lty = 1,lwd = 3,main = "Forecast of x3")